//
// Created by feng on 19-3-27.
// PS: Define some functions that can be used.
//
#pragma once

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

#include <cmath>
#include <cstdio>

#include "check_cuda.h"
#include "global_constant.h"
#include "vector_arithmetic.cuh"

namespace global_func
{
    // Rotate: The first step(1 of 2) to transform the local coordinate to world coordinate.

    __host__ __device__ inline auto local2world_rotate(const float3& d_local,
                                                       const float3& aligned_normal) -> float3
    {
        if (fabs(aligned_normal.x) < Epsilon && fabs(aligned_normal.z) < Epsilon)
        {
            return d_local;
        }

        // TODO::添加优化，判断n与(0,1,0)、(1,0,0)谁的距离更远，使用距离远的那个

        // printf("-----------------test local2world_rotate
        // begins---------------------------------\n"); printf(" d_local = %f %f %f  \n
        // aligned_normal = %f %f %f\n",
        //         d_local.x, d_local.y, d_local.z, aligned_normal.x, aligned_normal.y,
        //         aligned_normal.z);
        float3 u, n, v;
        n = aligned_normal;
        n = normalize(n);

        u = cross(make_float3(0.0f, 1.0f, 0.0f),
                  n); // !!!!!!!!!!!!!!!!!!!这里，如果aligned_normal就等于(0,1,0)该怎么办？？？
        u = normalize(u);
        v = cross(u, n);
        v = normalize(v);

        float3 d_world = make_float3(d_local.x * u.x + d_local.y * n.x + d_local.z * v.x,
                                     d_local.x * u.y + d_local.y * n.y + d_local.z * v.y,
                                     d_local.x * u.z + d_local.y * n.z + d_local.z * v.z);
        // printf("d_local = %f %f %f, n = %f %f %f, u = %f %f %f, v = %f %f %f,  d_world = %f %f %f
        // \n",
        //         d_local.x, d_local.y, d_local.z, n.x, n.y, n.z,u.x, u.y, u.z, v.x, v.y,
        //         v.z,d_world.x, d_world.y, d_world.z);
        // printf("-----------------test local2world_rotate
        // ends---------------------------------\n");

        return d_world;
    }

    // Translate: The second step(2 of 2) to transform the local coordinate to world coordinate.
    template <typename T>
    __host__ __device__ inline auto translate(const T& d_in, const T translate_vector) -> T
    {
        return d_in + translate_vector;
    }

    // Rotate matrix:
    // |cos		0		sin| |origin.x|     |cos * origin.x + sin * origin.z |
    // |0		1		0  | |origin.y|  =  |           origin.y             |
    // |-sin	0		cos| |origin.z|     |-sin * origin.x + cos * origin.z|
    __host__ __device__ inline auto rotateY(const float3& origin, const float3& old_dir,
                                            const float3& new_dir) -> float3
    {
        int dir = (cross(old_dir, new_dir).y > 0)
                      ? 1
                      : -1; // when parallel to (0,1,0), sin = 0, no effect
        float cos = dot(old_dir, new_dir);
        float sin = dir * sqrtf(1 - cos * cos);
        return make_float3(cos * origin.x + sin * origin.z, origin.y,
                           -sin * origin.x + cos * origin.z);
    }

    __host__ __device__ inline auto refractRay(const float3 dir1, const float3 normal,
                                               float refract_, float theta_1, float theta_2)
        -> float3
    {
        // 参数分别为：微定日镜中点坐标、入射光向量、入射光径向偏移
        // 计算折射入玻璃的向量         r2 = （r1/e + n1*(cosθ1/e – cosθ3))  *   H/ cosθ3
        float x = (dir1.x / refract_ + normal.x * (cos(theta_1) / refract_ - cos(theta_2)));
        float y = (dir1.y / refract_ + normal.y * (cos(theta_1) / refract_ - cos(theta_2)));
        float z = (dir1.z / refract_ + normal.z * (cos(theta_1) / refract_ - cos(theta_2)));
        return make_float3(x, y, z);
    }
    // Judge whether the ray intersects with the parallelogram (ractangle)
    /* Reference: Fast, Minimum Storage Ray Triangle Intersection
     * ray equation: orig + dir * t  (origin: orig, direction: dir)
     * A/B/C are three vertexes of the parallelogram and A is the angle subtend diagonal.
     * In rectangle, A = PI/2
     * orig + dir * t = (1-u-v)A + uB + vC
     *                   |t|
     * ==> |-D  B-A C-A| |u| = orig - A
     *                   |v|
     * Now use Cramer's Rule and the mixed product of vectors to solve the system :
     * Make E1 = B - A, E2 = C - A, T = orig - A
     * According to Cramer's Rule:
     *
     *       t =   _____1_____  |T E1 E2|
     *              |-D E1 E2|
     *
     *       U =   _____1_____  |-D T E2|
     *              |-D E1 E2|
     *
     *       V =   _____1_____  |-D E1 T|
     *              |-D E1 E2|
     *
     *  Combine these three solutions together:
     *      |t|                 | T E1 E2|
     *      |u| =  _____1_____  |-D  T E2|
     *      |v|     |-D E1 E2|  |-D E1  T|
     *  According to the mixed product of vectors: |a b c| = a x b · c = -a x c · b
     *
     *  it can be changed in this way:
     *      |t|                     |T x E1 · E2|
     *      |u| =   ______1______   |D x E2 · T |
     *      |v|     |D x E1 · E2|   |T x E1 · D |
     *  Make P = D x E2
     *       Q = T x E1
     *  Get the final formula:
     *      |t|                |Q · E2|
     *      |u| =   ___1____   |P · T |
     *      |v|     |P · E1|   |Q · D |
     */
    __host__ __device__ inline auto rayParallelogramIntersect(const float3& orig, const float3& dir,
                                                              const float3& A, const float3& B,
                                                              const float3& C, float& t, float& u,
                                                              float& v) -> bool
    {
        float3 AB = B - A;
        float3 AC = C - A;
        // vector p
        float3 p = cross(dir, AC);
        // determinant
        float det = dot(AB, p);

        // if determinant is close to 0, ray lies in plane of parallelogram
        if (fabsf(det) < Epsilon)
        {
            return false;
        }

        float3 T = orig - A;

        // Calculate u and make sure u <= 1
        u = dot(T, p) / det;
        if (u < 0 || u > 1)
        {
            return false;
        }

        // vector q
        float3 q = cross(T, AB);

        // Calculate v and make sure v <= 1
        v = dot(dir, q) / det;
        if (v < 0 || v > 1)
        {
            return false;
        }

        t = dot(AC, q) / det;
        return t >= Epsilon;
    }

    // copy from cpu to gpu
    template <typename T> inline void cpu2gpu(T*& d_out, T*& h_in, const size_t& size)
    {
        if (d_out == nullptr)
        {
            checkCudaErrors(cudaMalloc((void**)&d_out, sizeof(T) * size));
        }
        checkCudaErrors(cudaMemcpy(d_out, h_in, sizeof(T) * size, cudaMemcpyHostToDevice));
    }

    // copy from gpu to cpu
    template <typename T> inline void gpu2cpu(T*& h_out, T*& d_in, const size_t& size)
    {
        if (h_out == nullptr)
        {
            h_out = new T[size];
        }
        checkCudaErrors(cudaMemcpy(h_out, d_in, sizeof(T) * size, cudaMemcpyDeviceToHost));
    }

    // After fixing the number of threads
    __host__ __device__ inline auto setThreadBlocks(dim3& nBlocks, int nThreads, long long int size,
                                                    bool threadFixed) -> bool
    {
        if (size > MAX_ALL_THREADS)
        {
            printf("size %lld max %ld  There are too many threads to cope with. Please use less "
                   "threads.\n",
                   size, MAX_ALL_THREADS);
            return false;
        }
        // printf("size=%lld\n", size);
        int block_lastDim = (size + nThreads - 1) / nThreads;
        // printf("block_lastDim: %d\n" ,block_lastDim);

        if (block_lastDim < MAX_BLOCK_SINGLE_DIM)
        {
            nBlocks.x = block_lastDim;
            nBlocks.y = 1;
            nBlocks.z = 1;
            // printf("1nBlocks: ( %d, %d, %d )\n", nBlocks.x, nBlocks.y, nBlocks.z);
            return true;
        }

        block_lastDim = (block_lastDim + MAX_BLOCK_SINGLE_DIM - 1) / MAX_BLOCK_SINGLE_DIM;

        if (block_lastDim < MAX_BLOCK_SINGLE_DIM)
        {
            nBlocks.x = MAX_BLOCK_SINGLE_DIM;
            nBlocks.y = block_lastDim;
            nBlocks.z = 1;
            // printf("2nBlocks: ( %d, %d, %d )\n", nBlocks.x, nBlocks.y, nBlocks.z);
            return true;
        }

        nBlocks.x = MAX_BLOCK_SINGLE_DIM;
        nBlocks.y = MAX_BLOCK_SINGLE_DIM;
        nBlocks.z = (block_lastDim + MAX_BLOCK_SINGLE_DIM - 1) / MAX_BLOCK_SINGLE_DIM;
        // printf("3nBlocks: ( %d, %d, %d )\n", nBlocks.x, nBlocks.y, nBlocks.z);
        return true;
    }

    // fix the number of threads
    __host__ __device__ inline auto setThreadBlocks(dim3& nBlocks, int& nThreads,
                                                    const size_t& size) -> bool
    {
        nThreads = (MAX_THREADS <= size) ? MAX_THREADS : size;
        return setThreadBlocks(nBlocks, nThreads, size, true);
    }

    __host__ __device__ inline auto getThreadID() -> unsigned long long int
    {
        // unique block index inside a 3D block grid
        const unsigned long long int blockId = blockIdx.x                            // 1D
                                               + blockIdx.y * gridDim.x              // 2D
                                               + gridDim.x * gridDim.y * blockIdx.z; // 3D
        // global unique thread index, block dimension uses only x-coordinate
        const unsigned long long int threadId = blockId * blockDim.x + threadIdx.x;
        // printf("blockid=%lld  threadid=%lld\n",blockId,threadId);
        return threadId;
    }

    __host__ __device__ inline auto angle2xyz(float2 d_angles) -> float3
    {
        // printf("theta=%f\n",d_angles.x);
        return make_float3(sinf(d_angles.x) * cosf(d_angles.y), cosf(d_angles.x),
                           sinf(d_angles.x) * sinf(d_angles.y));
    }

    // Unroll and roll the index and address
    __host__ __device__ inline auto unroll_index(int3 index, int3 matrix_size) -> int
    {
        return index.x * matrix_size.y * matrix_size.z + index.y * matrix_size.z + index.z;
    }

    __host__ __device__ inline auto unroll_index(int2 index, int2 matrix_size) -> int
    {
        return index.x * matrix_size.y + index.y;
    }

} // namespace global_func
